home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Users Group Library 1996 July / C-C++ Users Group Library July 1996.iso / listings / v_08_05 / 8n05117a < prev    next >
Text File  |  1990-04-17  |  6KB  |  197 lines
  1. #include "math.h"
  2. #include "stdlib.h"
  3.  
  4. double adapt(x,y,n,res,err,lambda,numlambda,norm)
  5. /*
  6.  *  This routine adaptively fits a Hermite cubic to the data points
  7.  *  x[i],y[i] i = 0,1,...,n-1 according to the value of norm.
  8.  *          norm = 1  L1 (Robust) fit,
  9.  *               = 2  L2 (Least Squares) fit,
  10.  *               = 3  Minimax fit.
  11.  *  Output are the resulting coefficients in res[] and the errors
  12.  *  err[] at each point. The final fit is achieved using numlambda
  13.  *  knots lambda[].
  14.  *  The return value is the final value of the norm.
  15.  */
  16. double *x,*y,*res,*err,*lambda;
  17. int n,*numlambda,norm;
  18. {
  19. double m,capm,xm,*newlambda,*s,z,aic1,aic2;
  20. int j1,j2,i,k,l,newknots,numknots,flag=0;
  21. extern double (*faic)();
  22. /*  Allocate space */
  23. newlambda = (double*)calloc(2*n,sizeof(double));
  24. s = newlambda + n;
  25. /*  Fit two knots  */
  26. newlambda[0] = x[0]; newlambda[1] = x[n-1];
  27. z = Hermite (x,y,n,norm,newlambda,2,flag,res,err);
  28. aic1 = AIC(err,n,4,norm);
  29. /*  Fit three knots  */
  30. newlambda[2] = newlambda[1]; newlambda[1] = point(newlambda,res,1);
  31. z = Hermite(x,y,n,norm,newlambda,3,flag,res,err);
  32. aic2 = AIC(err,n,6,norm);
  33. numknots = newknots = 3;
  34. while (aic1>aic2)   /*  Test for minimum AIC */
  35.   {
  36.     aic1 = aic2;
  37.     for (i=0; i<numknots; i++) lambda[i] = newlambda[i];
  38.     numknots = newknots;
  39.     j2 = 0;
  40.     for (capm=0.0,i=0; i<n; i++) capm += err[i]*err[i];
  41.     capm /= (double)n;
  42.     /* 
  43.      *  Test each interval (lambda[k],lambda[k+1]) 
  44.      */
  45.     for (k=0; k<numknots-1; k++)
  46.       {
  47.         j1 = j2;
  48.         l = numpoints (x,n,lambda[k],lambda[k+1],&j2);
  49.         for (m=0.0,i=j1; i<j2+1; i++) m += err[i]*err[i];
  50.         m /= l;
  51.         if (m>capm)  /* Candidate point */
  52.           {
  53.             xm = point(lambda,res,k+1);
  54.             if (test_point(x,n,lambda,k) || k==0 || k-1==numknots)
  55.               {
  56.               /*  Add the knot to newlambda */
  57.                 for (i=k+1; i<newknots; i++)
  58.                 newlambda[i+1] = newlambda[i];
  59.                 newknots++;
  60.                 newlambda[k+1] = xm;
  61.               }
  62.           }
  63.       }
  64.     z = Hermite (x,y,n,norm,newlambda,newknots,flag,res,err);
  65.     aic2 = AIC(err,n,2*newknots,norm);
  66.   }
  67. /*  Fit is complete now optimize on minimum AIC */
  68. z = Hermite(x,y,n,norm,lambda,numknots,flag,res,err);
  69. for (i=1; i<numknots-2; i++) 
  70.   s[i] = log((lambda[i+1]-lambda[i])/(lambda[i]-lambda[i-1]));
  71. z = QN_BFGS_f(faic,s,numknots-3,err);
  72. return (z);
  73. }
  74.  
  75.  
  76.  
  77. double AIC (obs,n,p,flag)
  78. /*
  79.  *  This function computes the Akaike Information Criterion for 
  80.  *  the n observation obs[] with p paramaters. The distribution
  81.  *  assumed is given by flag, where
  82.  *        flag = 1  Cauchy,
  83.  *             = 2  Normal,
  84.  *             = 3  Minimax
  85.  */
  86. double *obs;
  87. int n,p,flag;
  88. {
  89. double xn,xp,sum,a,aic,xx,xbar; 
  90. int i;
  91. xn = (double)n; xp = (double)p;
  92. switch (flag)
  93.   {
  94.     case 1:
  95.       for (sum=0.0, i=0; i<n; i++) sum += obs[i]*obs[i];
  96.       a = sqrt(sum/xn);
  97.       for (aic=0.0, i=0; i<n; i++) aic += log(a*a + obs[i]*obs[i]);
  98.       aic -= (xn*log(a) - xp);
  99.       break;
  100.     case 2:
  101.       for (xbar=0.0,sum=0.0,i=0; i<n; i++)
  102.         { sum += obs[i]*obs[i]; xbar += obs[i]; }
  103.       xx = sum/xn - (xbar/xn)*(xbar/xn);
  104.       aic = xn*log(xx) + 2.0*xp;
  105.       break;
  106.     default:
  107.       for (sum=0.0,i=0; i<n; i++) sum += pow(obs[i],20.0);
  108.       aic = xn*log(sum)/20.0 + 2.0*xp;
  109.       break;
  110.   }
  111. return aic;
  112. }
  113.  
  114. double point(lambda,res,i)
  115. /*
  116.  *  This function computes the possible placement of an extra knot.
  117.  *  The present set of 'numknots' knots is given by lambda[] with a
  118.  *  present set of results res[]. The interval lambda[i-1],lambda[i]
  119.  *  is tested.
  120.  */
  121. double *lambda,*res;
  122. int i;
  123. {
  124. double a,b,c,dist,x[3],s,middle,xm,half,p;
  125. int k,k1,k2,k3,k4;
  126. k1 = 2*i - 2; k2 = k1 +1; k3 = k2 + 1; k4 = k3 +1;
  127. dist = lambda[i] - lambda[i-1];
  128. /*  
  129.  *  Compute derivative of cubic 
  130. */
  131. a = 6.0*res[k1] + 3.0*res[k2]*dist - 6.0*res[k3] + 3.0*res[k4]*dist;
  132. b = -6.0*(dist + 2.0*lambda[i-1])*res[k1];
  133. b -= 2.0*dist*(lambda[i-1] + 2.0*lambda[i])*res[k2];
  134. b += 6.0*(dist + 2.0*lambda[i-1])*res[k3];
  135. b -= 2.0*dist*(lambda[i] + 2.0*lambda[i-1])*res[k4];
  136. c = 6.0*lambda[i-1]*(dist + lambda[i-1])*res[k1];
  137. c += dist*lambda[i]*(lambda[i] + 2.0*lambda[i-1])*res[k2];
  138. c -= 6.0*lambda[i-1]*(dist + lambda[i-1])*res[k3];
  139. c += dist*lambda[i-1]*(lambda[i-1] + 2.0*lambda[i])*res[k4];
  140. /*
  141.  *   compute new knot
  142.  */
  143. x[0] = lambda[i-1]; x[1] = x[0]; 
  144. s = b*b - 4.0*a*c;
  145. if (s > 0.0)
  146.   {
  147.     s = sqrt(s); 
  148.     x[0] = (-b + s)/(2.0*a);
  149.     x[1] = c/x[0];
  150.   }
  151. x[2] = -b/(2.0*a);
  152. half = (lambda[i] - lambda[i-1])/2.0;
  153. p = middle = lambda[i-1] + half;
  154. for (k=0; k<3; k++)
  155.   if ((xm = fabs(middle - x[k])) < half) { half = xm; p = x[k]; }
  156. return p;
  157. }
  158.  
  159. numponts (x,n,xl,xr,k)
  160. /*
  161.  *  Compute the number of values from the data set x[n] in
  162.  *  the interval [xl,xr] starting with x[k]. Output the last
  163.  *  subscript in k.
  164.  */
  165. double *x,xl,xr;
  166. int n,*k;
  167. {
  168. int count = 0, j = *k;
  169. while (x[j] < xr)
  170.   if (x[j] > xl)
  171.     { count++; j++; }
  172. *k = j-1;
  173. return count;
  174. }
  175.  
  176. test_point (x,n,lambda,xm,k)
  177. /*
  178.  *  Routine to test whether the point xm is to be added to the 
  179.  *  new set of knots in the interval (lambda[k],lambda[k+1]).
  180.  */
  181. double *x,*lambda,xm;
  182. int n,k;
  183. {
  184. double x1,x2,x3;
  185. int insert,l1,l2,j=0;
  186. /*  Count points in the interval */
  187. l1 = numpoints (x,n,lambda[k],xm,&j);
  188. l2 = numpoints (x,n,xm,lambda[k+1],&j);
  189. x1 = xm - lambda[k]; x2 = lambda[k+1] - xm;
  190. x3 = ((x1<x2) ? x1 : x2)/(x1 + x2);
  191. insert = (x3 > 0.1 && ((x1<x2)? l1 : l2) > 10) ? 1 : 0;
  192. if (!insert)   /*  Exchange the knots  */
  193.   if (x1<x2) lambda[k] = xm;
  194.   else lambda[k+1] = xm;
  195. return insert;
  196. }
  197.